function calculateStatistics(Z,sites,triangles,gridStructure)
Z = Z(:,1:3);
nArcs = (size(Z,1)-1)/2; % number of arcs
N = 50; % granularity of lineary samplings along arc
time = 0:1/N:1;

totDist = 0; K = []; heights = [];
for i=1:nArcs
    ind = 2*(i-1)+1;
    disp(['Calculating statistics for arc = ' num2str(i) '/' num2str(nArcs)]);
    for j=1:length(time)-1
        t1 = time(j); t2 = time(j+1);
        
        % calculate distance
        F1 = ((1-t1)^2)*Z(ind,:) + 2*t1*(1-t1)*Z(ind+1,:) + (t1^2)*Z(ind+2,:);
        F2 = ((1-t2)^2)*Z(ind,:) + 2*t2*(1-t2)*Z(ind+1,:) + (t2^2)*Z(ind+2,:);
        totDist = totDist + norm(F1-F2); 
        
        % calculate curvature
        d = -2*(1-t1)*Z(ind,:) + 2*(1-2*t1)*Z(ind+1,:) + 2*t1*Z(ind+2,:);
        dd = 2*Z(ind,:) - 4*Z(ind+1,:) + 2*Z(ind+2,:);
        
        num = ((dd(3)*d(2)-dd(2)*d(3))^2 + (dd(1)*d(3)-dd(3)*d(1))^2 + (dd(2)*d(1)-dd(1)*d(2))^2)^(1/2);
        den = (d(1)^2 + d(2)^2 + d(3)^2)^(3/2);
        K = [K; num/den];
        
        % calculate safety
        triangleIndex = searchTriangle(F1,sites,triangles,1,gridStructure);
        if(triangleIndex==0)
            disp('Error finding Triangle. Location skipped')
        else
            nodeA = sites(triangles(triangleIndex,1),:);
            nodeB = sites(triangles(triangleIndex,2),:);
            nodeC = sites(triangles(triangleIndex,3),:);
            
            heights = [heights; checkSafety(F1, nodeA, nodeB, nodeC)];
        end
    end
end
disp(['Total Distance = ' num2str(totDist)]);

disp(['Maximum Curvature = ' num2str(max(K))]);
disp(['Average Curvature = ' num2str(mean(K))]);
disp(['Number of Arcs used = ' num2str(nArcs)]);

disp(['Minimum distance to terrain = ' num2str(min(heights))]);
disp(['Average distance to terrain = ' num2str(mean(heights))]);

stats = [totDist max(K) mean(K) nArcs min(heights) mean(heights)];

% save('../data/finalData/hw4.stats','stats','-ASCII', '-tabs')